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Abstract 

The problem of polynomial regression in which the usual monomial basis is re- 
placed by the Bernstein basis is considered. The coefficient matrix A of the overde- 
termined system to be solved in the least squares sense is then a rectangular 
Bernstein- Vandermonde matrix. 

In order to use the method based on the QR decomposition of A, the first stage 
consists of computing the bidiagonal decomposition of the coefficient matrix A. 

Starting from that bidiagonal decomposition, an algorithm for obtaining the QR 
decomposition of A is the applied. Finally, a triangular system is solved by using 
the bidiagonal decomposition of the R- factor of A. 

Some numerical experiments showing the behavior of this approach are included. 
AMS classification: 65F20; 65F35; 15A06; 15A23 
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1 Introduction 



Given {xj}i<.;<; + i pairwise distinct real points and {fi}i<i<i+i G R, let us 
consider a degree n polynomial 

P{x) = c + c±x + . . . + c n x n 
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for some n < I. Such a polynomial is a least squares fit to the data if it 
minimizes the sum of the squares of the deviations from the data, 

El^)-/,l 2 - 
i=i 

Computing the coefficients Cj of that polynomial P(x) is equivalent to solve, 
in the least squares sense, the overdetermined linear system Ac = f, where A 
is the rectangular (l + l) x (n + 1) Vandermonde matrix corresponding to the 
nodes {xj}i<i< J+ i. 

Taking into account that A has full rank n + 1, the problem has a unique 
solution given by the unique solution of the linear system 

A T Ac = A T f, 

the normal equations. 

Since A is usually an ill-conditioned matrix, it was early recognized that solv- 
ing the normal equations was not an adequate method. Golub [7], following 
previous ideas by Householder, suggested the use of the QR factorization of 
A, which involves the solution of a linear system with the triangular matrix 
R. 

Let us observe that, if A = QR with Q being an orthogonal matrix, then using 
the condition number in the spectral norm we have 

« 2 (i2) = K2(A), 

that is, R inherits the ill-conditioning of A while k 2 (A t A) = k 2 (A) 2 . 

In addition, as it was already observed by Golub in [8] (see also Section 20.1 
of [9]), although the use of the orthogonal transformation avoids some of the 
ill effects inherent in the use of normal equations, the value k 2 (A) 2 is still 
relevant to some extent. 

Consequently a good idea is to use, instead of the monomial basis, a polyno- 
mial basis which leads to a matrix A with smaller condition number than the 
Vandermonde matrix. 

It is frequently assumed that this happens when bases of orthogonal polyno- 
mials, such as the basis of Chebyshev polynomials, are considered. However, 
this fact is true when special sets of nodes are considered, but not in the case 
of general nodes. A basis which leads to a matrix A better conditioned than 
the Vandermonde matrix is the Bernstein basis of polynomials, a widely used 
basis in Computer Aided Geometric Designed due to the good properties that 
it possess (see, for instance, [2,10]). We illustrate these facts with Table 1, 
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where the condition numbers of Vandermonde, Chebyshev-Vandermonde and 
Bernstein- Vandermonde matrices are presented for the nodes considered in 
Example 5.1 and in Example 5.2. 



Example 


V 


TV 


BV 


5.1 


1.7e+12 


1.7e+12 


2.0e+05 


5.2 


2.5e+14 


4.0e+14 


5.3e+08 



Table 1 

Condition numbers of the Vandermonde V, Chebyshev-Vandermonde TV and 
Bernstein- Vandermonde BV matrices 



Let us observe that, without lost of generality, we can consider the nodes 
{ x i}i<i<i+i ordered and belonging to (0,1). So, we will solve the following 
problem: 

Let {xi}\<i<i + i G (0, 1) a set of points such that < x\ < . . . < xi+i < 1. Our 
aim is to compute a polynomial 

n 

P(x)= 5>&S n) (s) 

expressed in the Bernstein basis of the space H n (x) of the polynomials of 
degree less that or equal to n on the interval [0, 1] 

B n = [bf\x) = Q (1 - xf-V, j = 0, . . . , n}, 

such that P(x) minimizes the sum of the squares of the deviations from the 
data. 



This problem is equivalent to solve in the least squares sense the overdeter- 
mined linear system Ac = f where now 



A 



n-1 



(o)(l-^)" QxaCl-xa)- 1 



(»)^ 



l(s)( 1 - a; i+i) B 0^1+1(1 -xi+o"- 1 ••• Or + J 



is the (/ + 1) x (n + 1) Bernstein- Vandermonde matrix for the Bernstein basis 
i3 n and the nodes {xi}i<»<j+i, 



/ — (/l) /2, • • • , //+i) 



1.2) 
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is the data vector, and 

c = (ci,c 2 , . . . ,c n+1 ) T (1.3) 

is the vector containing the coefficients of the polynomial that we want to 
compute. 

The rest of the paper is organized as follows. Neville elimination and total 
positivity are considered in Section 2. In Section 3, the bidiagonal factorization 
of a rectangular Bernstein- Vandermonde matrix is presented. The algorithm 
for computing the regression polynomial in Bernstein basis is given in Section 
4. Finally, Section 5 is devoted to illustrate the accuracy of our algorithm by 
means of some numerical experiments. 



2 Basic results on Neville elimination and total positivity 

In this section we will briefly recall some basic results on Neville elimination 
and total positivity which we will apply in Section 3. Our notation follows the 
notation used in [4] and [5]. Given k, n G N (1 < k < n), Qk, n will denote the 
set of all increasing sequences of k positive integers less than or equal to n. 

Let A be an I x n real matrix. For k < I, m < n, and for any a G Qk,i and 
£ Qm,n, we will denote by the submatrix k x m of A containing the 

rows numbered by a and the columns numbered by (3. 

The fundamental tool for obtaining the results presented in this paper is the 
Neville elimination [4,5], a procedure that makes zeros in a matrix adding to 
a given row an appropriate multiple of the previous one. We will describe the 
Neville elimination for a matrix A = (o>i,j)i<i<i;i<j<n where I > n. The case in 
which I < n is analogous. 

Let A = (o>i,j)i<i<i;i<j<n be a matrix where / > n. The Neville elimination of 
A consists of n — 1 steps resulting in a sequence of matrices A := A\ — > A 2 — > 
. . . — > A n , where A t = (afj)i<i<i-i<j< n has zeros below its main diagonal in 
the t — 1 first columns. The matrix A t+1 is obtained from A t (t — 1, . . . , n) by 
using the following formula: 

afj , if i < t 

«S - (4V«ti>?^ , if *>t + landj>t + l (2.1) 
, otherwise. 

In this process the element 

Pi,j ■= <4j 1 < 3 < n , 3 < i < 1 



(*+i) ._ J 
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is called pivot of the Neville elimination of A. The process would break 
down if any of the pivots Pij (1 < j < n, j < % < I) is zero. In that case 
we can move the corresponding rows to the bottom and proceed with the new 
matrix, as described in [4]. The Neville elimination can be done without row 
exchanges if all the pivots are nonzero, as it will happen in our situation. The 
pivots pi t i are called diagonal pivots. If all the pivots pij are nonzero, then 
Pi,i = a?,iVi and, by Lemma 2.6 of [4] 



det A[i - j + 1, ... ... ,j] 

detA[i-j + l,...,i-l\l,...,j-l] 



Pij= ■ , v, ! : n1 Kj<n,j<i<l. (2.2) 



The element 

mj = -^- l<j<n, j<i<l (2.3) 

Pi-lJ 

is called multiplier of the Neville elimination of A. The matrix U := A n is 
upper triangular and has the diagonal pivots in its main diagonal. 

The complete Neville elimination of a matrix A consists on performing the 
Neville elimination of A for obtaining U and then continue with the Neville 
elimination of U T . The pivot (respectively, multiplier) of the complete 

Neville elimination of A is the pivot (respectively, multiplier) (j, i) of the 
Neville elimination of U T , if j >i. When no row exchanges are needed in the 
Neville elimination of A and U T , we say that the complete Neville elimination 
of A can be done without row and column exchanges, and in this case the 
multipliers of the complete Neville elimination of A are the multipliers of the 
Neville elimination of A if % > j and the multipliers of the Neville elimination 
of A T if j > i. 

A matrix is called totally positive (respectively, strictly totally positive) if all 
its minors are nonnegative (respectively, positive). The Neville elimination 
characterizes the strictly totally positive matrices as follows [4] : 

Theorem 2.1. A matrix is strictly totally positive if and only if its complete 
Neville elimination can be performed without row and column exchanges, 
the multipliers of the Neville elimination of A and A T are positive, and the 
diagonal pivots of the Neville elimination of A are positive. 

It is well known [2] that the Bernstein- Vandermonde matrix is a strictly totally 
positive matrix when the nodes satisfy < x\ < x 2 < . . . < xi +1 < 1, but this 
result will also be shown to be a consequence of our Theorem 3.2. 
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3 Bidiagonal factorization of A 



In this section we consider the bidiagonal factorization of the Bernstein- 
Vandermonde matrix A of (1.1). 

Let us observe that when I = n this matrix A is the coefficient matrix of 
the linear system associated with a Lagrange interpolation problem in the 
Bernstein basis B n whose interpolation nodes are {xi : i — 1, . . . , n + 1}. A 
fast and accurate algorithm for solving this linear system, and therefore the 
corresponding Lagrange interpolation problem in the Bernstein basis can be 
found in [13]. A good introduction to the interpolation theory can be seen in 



The following two results will be the key to construct our algorithm. 

Proposition 3.1. (See [13]) Let A be the square Bernstein- Vandermonde 
matrix of order n + 1 for the Bernstein basis B n and the nodes x±, x 2 , ■ ■ ■ , x n+ \. 
We have: 



Theorem 3.2. Let A = (aij)i<i<;+i ; i<j< n +i be a Bernstein- Vandermonde 
matrix for the Bernstein basis B n whose nodes satisfy < x 1 < x 2 < . . . < 
x\ < xi + \ < 1. Then A admits a factorization in the form 



where Gj are (n + 1) x {n + 1) upper triangular bidiagonal matrices (j = 
1, . . . ,n), Fi are (/ + 1) x (/ + 1) lower triangular bidiagonal matrices (i = 
1, . . . , /), and D is a (I + 1) x (n + 1) diagonal matrix. 

Proof. The matrix A is strictly totally positive (see [2]) and therefore, by 
Theorem 2.1, the complete Neville elimination of A can be performed without 
row and column exchanges providing the following factorization of A (see [6]): 



[3]. 




A = F l F l . 1 ---F 1 DG 1 ---G n - 1 G.. 



(3.1) 



A = F l F l _ 1 ---F 1 DG 1 ---G n - 1 G. l 
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where Fj (1 < i < I) are (I + 1) x (/ + 1) bidiagonal matrices of the form 



Fi = 



1 

1 



1 

m i+2 ,2 1 



V 



mii-i 1 



(3.2) 



Gj (1 < « < n) are (n + 1) x (n + 1) bidiagonal matrices of the form 



G 



1 

1 



1 

rn i+lj i 1 

m+2,2 l 



V 



m r , 



(3.3) 
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and D is the (/ + 1) x (n + 1) diagonal matrix whose ith (1 < % < n + 1) 



diagonal entry is the diagonal pivot pij = afj of the Neville elimination of A: 



D — (rfij)l<i<i+l;l<j<n+l — diag{pi ; i,p 2 ,2, • • • ,Pn+l,n+l}- 



(3.4) 



Taking into account that the minors of A with j initial consecutive columns 
and j consecutive rows starting with row i are 



detA[i,...,i + j - 



(1 _ Xi )»-j+l(l - x i+1 ) n -* +1 • • • (1 - Xi+j-J"-* 1 1] (*h ~ X k ), 

i<k<h<i+j-l 



a result that follows from the properties of the determinants and Proposition 
3.1, and that are the multipliers of the Neville elimination of A, we obtain 
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that . ^ 

Pi,j (1 x i) 3 (1 x i—j) rifc=l( a 'i x i— k) /q r\ 

(i-^-i)^' +2 ni =2 (^-i-^) ' 

where j = 1, . . . , n + 1 and i — j + 1, . . . , I + 1. 

As for the minors of A T with j initial consecutive columns and j consecutive 
rows starting with row i, they are: 

det^[v.. I+ ,-ni,...,^( j : i )Q-( j+ ;_ 2 ),r'^-4- 1 
(i-x 1 r-'-i+Hi-x 2 r— i+2 ---(i'x i r-'- i+2 n (*«.-**)• 

l<k<h<j 

This expression also follows from the properties of the determinants and 
Proposition 3.1. Since the entries mjj are the multipliers of the Neville elimi- 
nation of A T , using the previous expression for the minors of A T with initial 
consecutive columns and consecutive rows, it is obtained that 

_ (n — i + 2) • Xj . . 

mi -i = r — V j = 1, . . . , n; i = j + 1, . . . , n + 1. (3.6) 

Finally, the ith diagonal element of D 



Pa — - — : , i — l,...,n + l (3.7) 

ir fc =i(i-a*) 

is obtained by using the expression for the minors of A with initial consecutive 
columns and initial consecutive rows. □ 

Moreover, by using the same arguments of [14] it can be seen that this factor- 
ization is unique among factorizations of this type, that is to say, factorizations 
in which the matrices involved have the properties shown by formulae (3.2), 
(3.3) and (3.4). 



Remark 3.3. The formulae obtained in the proof of Theorem 3.3 for the 
minors of A with j initial consecutive columns and j consecutive rows, and for 
the minors of A T with j initial consecutive columns and j consecutive rows 
show that they are not zero and so, the complete Neville elimination of A 
can be performed without row and column exchanges. Looking at equations 
(3.5)-(3.7) is easily seen that m it j, and p iti are positive. Therefore, taking 
into account Theorem 2.1, this confirms that the matrix A is strictly totally 
positive. 



Remark 3.4. In the square case, the matrices Fi (i = 1, . . . , I) and the ma- 
trices Gj (j = 1, . . . ,n) are not the same bidiagonal matrices that appear in 
the bidiagonal factorization of A^ 1 presented in [13], nor their inverses. The 
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multipliers of the Neville elimination of A and A T give us the bidiagonal fac- 
torization of A and A^ 1 , but obtaining the bidiagonal factorization of A from 
the bidiagonal factorization of A^ 1 (or vice versa) is not straightforward [6]. 
The structure of the bidiagonal matrices that appear in both factorizations 
is not preserved by the inversion, that is, in general, F^ 1 (i = 1, . . . ,/) and 
Gj 1 (j = l,...,n) are not bidiagonal matrices. See [6] for a more detailed 
explanation. 



4 The algorithm 

In this section we present an accurate and efficient algorithm for solving the 
problem of polynomial regression in Bernstein basis we have presented in Sec- 
tion 1. As we introduced there, our algorithm is based on the solution of the 
least squares problem min c \\ Ac — f ||, where A, f and c are given by (1.1), 
(1.2) and (1.3), respectively. Taking into account that A is a strictly totally 
positive matrix, it is full rank, and the method based on the QR decomposition 
is the most adequate for solving this least squares problem [1]. 

The following result (see Section 1.3.1 in [1]) will be essential in the construc- 
tion of our algorithm. 

Theorem 4.1. Let Ac = f a linear system where A G R( z + 1 ) x ( n + 1 ) ) / > n , 
c G R n+1 and / G R m . Assume that rank(A) — n + 1, and let the QR 
decomposition of A given by 



where Q G R('+ 1 ) x ('+ 1 ) is an orthogonal matrix and R G R,( n + 1 ) x ( n + 1 ) [ s an 
upper triangular matrix with nonnegative diagonal entries. 

The solution of the least squares problem min c || Ac — f || 2 is obtained from 



where d\ G R n+1 , d 2 G R' n and r — f — Ac. In particular || r H2HI ^2 lb- 

An accurate and efficient algorithm for computing the QR decomposition of a 
strictly totally positive matrix A is presented in [12]. This algorithm is called 
TNQR and can be obtained from [11]. Given the bidiagonal factorization of A, 





Rc = d 
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TNQR computes the matrix Q and the bidiagonal factorization of the matrix 
R. Let us point out here that if A is strictly totally positive, then R is strictly 
totally positive. TNQR is based on Givens rotations, has a computational cost of 
0(l 2 n) arithmetic operations if the matrix Q is required, and its high relative 
accuracy comes from the avoidance of subtractive cancellation. 

A fast and accurate algorithm for computing the bidiagonal factorization of 
the rectangular Bernstein- Vandermonde matrix that appears in our problem 
of polynomial regression in the Bernstein basis can be developed by using the 
expressions (3.5), (3.6) and (3.7) for the computation of the multipliers rriij 
and m it j, and the diagonal pivots p iti of its Neville elimination. The algorithm 
is an extension to the rectangular case of the one presented in [13] for the 
square Bernstein- Vandermonde matrices. Given the nodes {xi}±<i<i + i G (0, 1) 
and the degree n of the Bernstein basis, it returns a matrix M £ R('+ 1 ) x ( n+1 ) 
such that 



The algorithm, that we call it TNBDBV, has a computational cost of 0{ln) arith- 
metic operations, and high relative accuracy because it only involves arith- 
metic operations that avoid subtractive cancellation (see [13] for the details). 
The implementation in Matlab of the algorithm in the square case can be 
taken from [11]. 

In this way, the algorithm for solving the least squares problem min c \\ Ac—f || 
corresponding to our polynomial regression problem will be: 

INPUT: The nodes {xi}i<i<i+i G (0, 1), the data vector / and the degree n of 
the Bernstein basis. 

OUTPUT: A vector c = {cj)-i<i< n+1 containing the coefficients of the polyno- 
mial P(x) in the Bernstein basis B n and the minimum residual r. 

- Step 1: Computation of the bidiagonal factorization of A by means of 
TNBDBV. 

- Step 2: Given the matrix M obtained in Step 1, computation of the QR 
decomposition of A by using TNQR. 

- Step 3: Computation of 




Pi,i 



m,j j = 1, • • • , n + 1; % = j + 1, ...,/ + 1, 
% % = 1, . . . , n\ j = % + 1, . . . , n + 1. 



i = 



l,...,n + l, 




- Step 4'- Solution of the upper triangular system Rc = d\. 
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Step 5: Computation of 




Step 3 and Step 5 are carried out by using the standard matrix multiplication 
command of MATLAB. As for Step 4, it is done by means of the algorithm 
TNSolve of P. Koev [11]. Given the bidiagonal factorization of a totally positive 
matrix A, TNSolve solves a linear system whose coefficient matrix is A by using 
backward substitution. 

Let us observe that A is not constructed, although we are also computing the 
residual r — f — Ac. 



5 Numerical experiments and final remarks 



Two numerical experiments illustrating the good properties of our algorithm 
are reported in this section. We solve the least squares problem min c \\ Ac—f || 
corresponding to the computation of the regression polynomial in exact arith- 
metic by means of the command least sqrs of Maple 10 and we denote this 
solution by c e . We also compute the minimum residual r e in exact arithmetic 
by using Maple 10. We use c e and r e for comparing the accuracy of the results 
obtained in MATLAB by means of: 

(1) The algorithm presented in Section 4. 

(2) The command A\f of Matlab. 

The relative errors obtained when using the approaches (1) and (2) for comput- 
ing the coefficients of the regression polynomial in the experiments described 
in this section (eci and ec2, respectively) are included in the first and in the 
third column of Table 2. The relative errors corresponding to the computation 
of the minimum residual by using the approaches (1) and (2) (eri and er 2 , 
respectively) are presented in the second and in the fourth column of Table 2. 

We compute the relative error of a solution c of the least squares problem 
min c || Ac — f || by means of the formula 

|| C - C e || 2 

ec = — 7 r . 

II C e 1 1 2 

The relative error of a minimum residual r is computed by means of 



Example 5.1. Let B15 the Bernstein basis of the space of polynomials with 
degree less than or equal to 15 in [0, 1]. We will compute the polynomial 

P{x) = Y,c 3 bf\x) 

3=0 

that minimizes 
where 

{^i}l<i<21 = <^ — \ , 
I ZZ J l<j<21 

and 

/ = (3, 4, 0, -2, 5, 0, 1, 9, -3, 7, -1, 0, 2, 2, -4, -2, 3, 8, -6, 4, if. 

Let us observe that, the condition number of the Bernstein- Vandermonde ma- 
trix A of the least squares problem corresponding to the regression polynomial 
we are interested in computing is k 2 (A) = 2.0e + 05. 

The following example shows how the algorithm we have presented in this 
paper keeps the accuracy when the condition number of the Bernstein- Vander- 
monde matrix involved in the regression problem increases, while the accuracy 
of the general approach (2) which does not exploit the structure of this matrix 
goes down. 

Example 5.2. We consider a regression problem such that the Bernstein 
basis £>i5 and the data vector / are the same as in Example 5.1. The points 
{xi}i<i<2i are now: 

f 1 1 1 1 1 1 1 1 1 1 1 23 21 19 17 15 13 11 9 7 51 

1 22' 20' 18' 16' 14' 12' 10' 8' 6' 4' 2' 42' 38' 34' 30' 26' 22' 18' 14' 10' 6 J" 

The condition number of the Bernstein- Vandermonde matrix A involved in 
this experiment is k 2 (A) = 5.3e + 08. 



Example 


eci 


eri 


ec 2 


er 2 


5.1 


1.4e-15 


1.3e-15 


9.0e-12 


5.2e-12 


5.2 


2.0e-15 


2.3e-15 


1.9e-09 


1.4e-08 



Table 2 

Relative errors in Example 5.1 and Example 5.2 

Remark 5.3. The accuracy of our algorithm is obtained by exploiting the 
structure of the Bernstein- Vandermonde matrix. Every step of our algorithm, 
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except the ones in which the standard matrix multiplication command of 
Matlab is used, are developed with high relative accuracy because only arith- 
metic operations that avoid subtractive cancellation are involved [12,13]. 

Remark 5.4. Our algorithm has the same computational cost (0(l 2 n) arith- 
metic operations) as the conventional algorithms that solve the least squares 
problem by means of the QR decomposition ignoring the structure of the 
matrix, when Q is explicitly required (see Section 2.4.1 of [1]). 
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